Collision statistics for random flights with anisotropic scattering and absorption 
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For a broad class of random walks with anisotropic scattering kernel and absorption, we derive 
explicit formulas that allow expressing the moments of the collision number nv performed in a 
volume y as a function of the particle equilibrium distribution. Our results apply to arbitrary 
domains V and boundary conditions, and allow assessing the hitting statistics for systems where 
the typical displacements are comparable to the domain size, so that the diffusion limit is possibly 
not attained. An example is discussed for one-dimensional (Id) random flights with exponential 
displacements, where analytical calculations can be carried out. 



I. INTRODUCTION 

The dynamics of complex physical systems is often de- 
scribed in terms of 'particles' undergoing random dis- 
placements, resulting either from the intrinsic stochas- 
tic nature of the underlying process, or from uncer- 
tainty [H, 0]. Widespread examples arise in radiation 
transport, research strategies, biology andpercolation 
through porous media, only to name a few [3|-|6|. In this 
context, quantifying the residence time that the walkers 
spend inside a given domain is a key issue that has 
motivated a considerable research effort 

When the domain size is large as compared to the aver- 
age displacements, the walker dynamics is usually mod- 
elled by either regular Brownian motion for homogenous 
media, or anomalous diffusion for heterogeneous, scale- 
invariant media P, S, i, For Markovian transport 
processes, a systematic approach to assessing the resi- 
dence time distribution exists, via the so-called Feynman- 
Kac formalism p"8l - [2^ . Yet, full knowledge of the resi- 
dence time distribution is an awkward task, and is achiev- 
able only in a limited number of cases [lH, , so that 
one has often to be content with the first few moments of 
the residence time 0, Bin, HI. A further difficulty arises 
when the walker typically undergoes a limited number of 
collisions before leaving the explored domain, and the dif- 
fusion limit is possibly not attained. This is often the case 
in gas dynamics, neutronics and radiative transfer, elec- 
tronics, and biology 5, 26r.29i] . In all such systems, the 
stochastic path can be thought of as a series of straight- 
line flights, separated by random collisions, and a natural 
variable for describing the walker evolution is therefore 
the number of collisions ny within the observed volume. 
Application of the diffusion approximation to the count- 
ing statistics, which amounts to assuming a large number 
of collisions in V, might lead to inaccurate results [13] ■ 

In a previous work, we have addressed the issue of 
characterizing the moments {riy) for arbitrary geome- 
tries and boundary conditions, subject to the condition 
that both the source and scattering are isotropic [Slj . 
Here we extend those results by relaxing the isotropy hy- 



pothesis. We also distinguish the case where events are 
counted before or after each collision. We derive explicit 
formulas for the moments (n™) by building on survival 
probabilities, and relate the collision statistics to the par- 
ticle equilibrium distribution. Knowledge of higher order 
moments allows estimating the uncertainty on the aver- 
age, as well as reconstructing the full distribution of the 
collision number. We exemplify our findings by examin- 
ing a case of one-dimensional (Id) transport, the so called 
rod model with exponentially distributed displacements, 
where analytical calculations can be carried out. 

This paper is structured as follows: in Sec.|lT]we briefly 
recall the basic properties of random flights performing 
anisotropic scattering and absorption. In Sec. IIIII we de- 
rive the moments (riy) by a direct contruction based on 
survival probabilities. Then, in Sec. IIVI we discuss the 
asymptotic results that are recovered in the large ny 
limit. Finally, an example is worked out in Sec. |Vl and 
conclusions are drawn in Sec. IVII 



II. 



RANDOM FLIGHTS: TRANSPORT 
KERNELS AND EQUILIBRIUM 
DISTRIBUTIONS 



Consider the random walk of a particle emitted at ve- 
locity Vq from a point-source S located at Tq. At each 
collision, the particle can be either scattered, with prob- 
ability p(r, v), or absorbed (in which case the trajectory 
terminates). 

We introduce the quantity T(r|r',v'), namely, the 
probability density of performing a displacement from 
r' to r (in the direction of v'), between any two colli- 
sions. Then, Jy r(r|r', v')dr represents the average num- 
ber of next collisions in a volume V per particle emitted 



at r', with velocity v' 



Analogously, we intro- 



duce C(v|v', r), namely, the conditional probability den- 
sity of changing velocity from v' to v, given a scatter- 
ing event at r. Usually, we can factorize C(v|v',r) = 
p(r, v')c(v|v', r), where c(v|v', r) denotes the normalized 
scattering kernel. Then, C(v|v', r')dv represents the 
average number of particles leaving a collision per inci- 
dent particle entering the collision at r', with velocity 



V 



[12, il]. It follows that 
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K{r, v|r', v') = r(r|r', v)C(v|v', r') 



(1) 
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FIG. 1: (Color online) Average collision number (n{/)p(a;o) 
for a sphere with radius R — 1 and leakage boundary condi- 
tions. The source is isotropic, andp = 1. Solid lines: Eq. (|45|l . 
Symbols: Monte Carlo simulations. The persistence parame- 
ter is a = (dark blue stars), a — 0.25 (green circles), a = 0.5 
(red squares), a = 0.75 (light blue dots), and a = 1 (violet 
diamonds) . 



FIG. 2: (Color online) Second moment of the collision number 
{nY)p{xo) for a sphere with radius R = 1 and leakage bound- 
ary conditions. The source is isotropic, andp — 1. Solid lines: 
Eq. (|47p . Symbols: Monte Carlo simulations. The persistence 
parameter is a = (dark blue stars), a = 0.25 (green circles), 
a = 0.5 (red squares), a = 0.75 (light blue dots), and a — 1 
(violet diamonds). 



represents the density of particles entering the (n-|- l)-th 
collision with coordinates r, v, having entered the n-th 
with coordinates r',v'. Inversing the order of the dis- 
placement and collision kernels, the quantity 

L(r, v|r', v') = C(v|v', r)r(r|r', v') (2) 

represents the density of particles leaving the 7i-th colli- 
sion with coordinates r, v, having left the (n— l)-th with 
coordinates r', v'. 

We can now define the displacement operator 

T[/](r,v')= J T(r|r',v')/(r',v')dr', (3) 

and the collision operator 

C[/](r,v)= J C{v\v',r)f{r,V)dV (4) 

for any sufficiently well-behaved /. Furthermore, we in- 
troduce the transport operators 

K[/](r,v) = Jdr' J dv'A-(r,v|r',v')/(r',v'), (5) 

and 

L[/](r,v) =Jdv' J dr'i(r,v|r',v')/(r',v'). (6) 

To simplify notation, we denote in the following z — 
{r, v} the coordinates of the walker in the phase space. 
We introduce then the incident propagator ^(z,r7,|zo), 
i.e., the probability density of finding a particle entering 



the n-th collision with coordinates z, starting from zq, 
and the outgoing propagator x(z, njzo), i.e., the probabil- 
ity density of finding a particle exiting the n-th collision 
with coordinates z, starting from zq. These quantities 
are related by 

x(z,n|zo) = C«'(z,n|zo) (7) 

and 

*(z,n|zo) = Tx(z,n- l|zo), (8) 
with n > 1, and 0|zo) = S. It follows also 

«'(z,n+l|zo) =K«'(z,n|zo) (9) 

and 

x(z,n|zo) = Lx(z,"- - l|zo)- (10) 

In other words, knowledge of the system state z at n is 
sufficient to determine the state at n + 1. From Eqs. © 
and (jlOp. by recursion we have 

^(z,n|zo) =K"-^T[5] (11) 

and 

x(z,n|zo) = L"[5], (12) 

where A" is a n-fold iterated operator. These relations 
show that the particle dynamics is entirely defined in 
terms of the two kernels C and T. 
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FIG. 3: (Color online) Average collision number (n^)p(n, a;o) 
for a sphere with radius R = 1 and leakage boundary condi- 
tions. The source is isotropic (xq = 0), and p = 1. Dashed 
lines: Eq. (|45p . Symbols: Monte Carlo simulations. The per- 
sistence parameter is a = (dark blue stars), a = 0.25 (green 
circles), a = 0.5 (red squares), a — 0.75 (light blue dots), and 
a — 1 (violet diamonds). 
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FIG. 4: (Color online) Second moment of the collision num- 
ber {riy) p{n, xo) for a sphere with radius R — 1 and leakage 
boundary conditions. The source is isotropic {xo = 0), and 
p — 1. Dashed lines: Eq. (|47p . Symbols: Monte Carlo simu- 
lations. The persistence parameter is a = (dark blue stars), 
a = 0.25 (green circles), a — 0.5 (red squares), a = 0.75 
(light blue dots), and a — 1 (violet diamonds). 



We introduce now the incident and outgoing collision 
densities, respectively. 



N 



*(z|zo) = Jim V4'(z,n|zo), 

n=l 
N 

x(z|zo) = lim Vx(z,n|zo), 



(13) 



n=0 



which can be interpreted as the particle stationary dis- 
tributions [3^ H^l • We can associate to the collision den- 
sities their respective operators, namely, 

vl>[/](z) = I vl/(z|z')/(z')rfz', 

X[/](z) = J x(z|z')/(z')dz'. (14) 

In particular, = 5'(z|zo), and x[S] = x(z|zo). 

Now, by making use of the formal Neumann series (see 
Eq. !(M\ ). from Eqs. dTT]) and ([H]) we have then 



T 



I-K 



and 



X[/](z) 



I 



(15) 



(16) 



Finally, it follows that the incident collision densities sat- 
isfies the stationary integral transport equation 



*(z|zo) =K*(z|zo)+T5 



(17) 



whereas the outgoing collision density satisfies 

x(z|zo) = Lx(z|zo)+5. (18) 
The solutions \E'(z|zo) and x(z|zo) are related by 



x(z|zo) = C*(z|zo)-l-5. 



(19) 



Observe that x(z|zo) obeys an integral equation whose 
source term is the physical source S, whereas the source 
term in the equation for ^'(z|zo) is the so-called first- 
collision source TS, i.e., the density of particles entering 
the first collision. For reasons that will be clear later, 
it is expedient to introduce the function (p(z\zi), being 
the solution of (l — K)'/3(z|zi) = S' , for a point-source 
consisting in a particle entering the first collision at Zi. 
Then, the incident collision density can be expressed by 
the convolution 

«'(z|zo) = / (^(z|zi)r(zi|zo)dzi. (20) 



III. COLLISION STATISTICS 

Suppose that the trajectories of the random flights de- 
scribed above are observed until the walker either disap- 
pears by leaving an external boundary, or is absorbed. 

A. Scattering and absorption events 

We first assume that each event in a given volume V 
is detected when the particle enters a collision (in other 
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FIG. 5: (Color online) The distribution P(yv) for a;o = 
as a function of the rescaled variable yv = nv / {rjR^). Left: 
71 = 30 and a = 0.25 (blue dots), a = 0.5 (green squares), and 
a = 0.75 (red diamonds). Right: a — 0.5 and R — 30 (blue 
dots), R = 40 (green squares), and -R = 50 (red diamonds). 



words, we do not discriminate scattering and absorption 
events). The quantity qqi{n\zQ) — Jy dz"^(z,n\zQ) rep- 
resents the survival probability, i.e., the probability for 
the particle to be in V up to entering the n-th collision. 
From the Markovian nature of the process z [sj , it fol- 
lows that the probability of detecting ny collision events 
in the volume V is 

P{nv\zo) ^ q^{nv\zo) - q^{nv + l\zo). (21) 

The moments are given by 

«)p(zo) - <Pinv\^o) (22) 

for m > 1, and depend on the boundary conditions on 
dV, which affect the functional form of the propaga- 
tor [35^. Setting boundary conditions at infinity corre- 
sponds to defining a fictitious ('transparent') volume V, 
where particles can indefinitely cross dV back and forth. 
On the contrary, the use of leakage boundary conditions, 
i.e., when the particle is lost upon crossing the boundary, 
leads to the formulation of first-passage problems 0, 

Normalization implies (ny)p(zo) = 1, and direct cal- 
culation from Eqs. and (PT|) yields 



(ny)p(zo) 



dz^[S] 



dz*(z|zo). (23) 



Observe that the integral of the collision density over a 
volume V gives the mean number of collisions within that 
domain, hence the name given to 4'(z|zo). 

Higher order moments of ny follow from Eqs. (jlip 



and (|22p . and read 

+ 00 



V. 



p(zo) = V < / dzK"--i(l-K)T[5], (24) 
Jv 



nv— 1 



for m > 1. By resorting to the identity ()A2p . we obtain 
then 

«)p(zo) = V klS^.k / rfz* (C*)'"' [S], (25) 
k=i -'^ 

where Sm,k are the Stirling numbers of second kind (see 
Eq. (jA3p ). We can further introduce the factorial mo- 



ments (n(7'^)(zo), where x^*^^ = x{x + l)...{x + k — 1) 
is the rising factorial [s^. The factorial moments are 
related to the moments by 



(rn) 



(26) 



A:=0 



Sm,k being the Stirling numbers of first kind (see 
Eq. (jASp ). From the operator identity (|A4p . we have 
then 

(n(7'))p(zo) - ml j^dz (^-^^^[S], (27) 
which holds for m > 1. Setting now 

T \™ 

[5'] (28) 



(4"))^(zi) = m! I^dz(^^ 



K 



yields the recursion property for the factorial moments 
(nlr^)i.(zi) = m / dz^(z|zi)(n(r-'Vp(z), (29) 



for rn > 1, starting from {riy )'p{zi) — 1. Then, from 
Eq. ((201) we have 

(zo) = J {n'(J'%{z,)T{z,\zo)dz,. (30) 



{ny )P 



B. Scattering events 

Suppose now that events in V are detected at the exit 
of each collision (in other words, only scattering col- 
lisions and the source S are recorded). The quantity 
(7x(n|zo) = dzx(z, n|zo), n > 0, represents the proba- 
bility for the particle to be in V after having undergone 
the n-th collision. From the same argument as above, 
the probability of detecting ny scattering events in the 
volume V is 

Q{nv\zQ) = q^{nv - l|zo) - ^^("vlzo)- (31) 

The moments are then 

«)q(zo) = <Q("v|zo) (32) 



5 



for m > 1, with (ny)Q (zq) = 1 from normalization. Di- 
rect calculation from Eqs. (15^ and (pij) yields 

(n^)Q(zo) = / dzx[S] = / dzx(z|zo). (33) 
Jv Jv 

Higher order moments of ny follow from Eq. , and 
read 

+ 00 „ 

«)q(zo)- V n'^ / dzL"--i(l-L)[5], (34) 
„ _i Jv 



which holds for m > 1. By resorting to the identity (|A2j), 
we obtain 

(n^)Q(zo) = Vfc!5™.fc / dzX(LX)''"' [5]. (35) 
Jv 



As done before, we can further introduce the factorial 
moments (n^^)(zo)Q. From the identity (|A4| . for m > 1 
we have then 



(nlr%(zo)=m! J^dz(^j^J 



[S], (36) 



which finally yields the recursion property for the facto- 
rial moments 



Equation (|39|) implies that at large n any anisotropic 
walk will behave like an isotropic walk (with a linear 
spread), provided that we rescale n by a factor 77 = 
1/ [1 -I- 27/i/(l — /i)] that depends on the specific features 
of the jump length distribution. For a sharply forward- 
peaked walk {fJ. 1), the number of collisions needed for 
attaining isotropy becomes very large as compared to an 
isotropic walk, for a given (^^). When n is sufficiently 
large, we expect the same scaling to carry over to the 
moments, i.e.. 



«)p(zo)o 



77'"«)p(zo),,o. 



(40) 



The so-called diffusion limit is reached when (^^) <C 
and the average flight time t = {£)/v is vanishing small, 
while preserving a constant ratio D — (£^)/r, namely 
the diffusion coefficient. Under these conditions, the 
collision number in V diverges, whereas the quantity 
ty — X]r=i converges to the residence time in the 
volume. Actually, ty should take into account also ad- 
ditional terms due to boundary conditions. However, as 
T — 0, the trajectory will almost surely have a turning 
point touching the boundary, so that corrections can be 
safely neglected. By following the arguments in [3l|, the 
moments (i™)(ro) of the residence time will be given by 
the celebrated Kac formula for the isotropic Brownian 
motion [HI, HI], up to a scaling factor 77™. 



("'y"%(zo) = "^ / dzx(z|zo)(n^" 'Oq(z). (37) 
Jv 

for m > 1. Similar results for the factorial moments 
appear in (37j (and references therein) under the name 
of discrete Feynman moment formula. 



IV. DIFFUSION LIMIT 

Suppose that the walker evolves in a medium without 
boundaries, starting from an isotropic source. Assume 
for the sake of simplicity that there is no absorption, and 
the speed v — |v| is constant. The spread at the n- 
th collision {r'^){n) = J dr J dvjr — rop^(r, v, Tijro, vo) 
reads 



{r')in) = n{f) + 2n(£)2-^ - 2{tf ^^^ (38) 
1 — /i 



(1-M)^ 



where I is the inter-collision length, and < /i < f is 
the average (polar) scattering angle [s^. Observe that 
Eq. (|38p does not explicitly depend on the dimension d 
of the embedding space. When the typical displacements 
are much smaller than the domain size, say ^ i?^, 
n is large, and the last term in Eq. ([38]) can be dropped. 
If the two first moments of the flight length distribution 
are finite, we can set 7 = {tf' j ij/?) > 0, and rewrite 



V. THE ROD MODEL 

The approach presented in Sec. |llT] allows explicitly 
evaluating the moments (ny)(ro, Vq). When the equilib- 
rium distribution is known, this amounts to solving the 
convolution integrals in Eqs. and ([55]) . However, 
analytical expressions for y (r, v|rn, Vp) can be obtained 
only in a few cases [3l|, HI, [H, Eol , so that one must gen- 
erally resort to numerical integration. In this Section we 
exemplify the moments formulas above for a well-known 
and long-studied system where calculations can be car- 
ried out analytically. 



A. Exponential flights 

When the scattering centers are spatially uniform, i.e., 
when the traversed medium is homogeneous, the inter- 
collision lengths are exponentially distributed. Exponen- 
tial flights describe for instance the displacements of neu- 
tral particles (neutrons, photons) in matter, the motion 
of electrons in semiconductors, the migration of biologi- 
cal species (often called velocity jump proc ess), an d are 
widely used in gas dynamics (Lorentz gas) |4l| - |47| . The 
displacement kernel for exponential flights reads 

r(r|r',v) = Et(r', v)e- ^*('"'+'*'^'^)'*^ (41) 



(39) 



when CO = v/v is parallel to r — r' , oriented as v [32, [33| . 
The quantity Sj(r,v) is the total cross section, which is 



6 



proportional to the probability of particle-medium inter- 
action along a straight line (it carries units of the inverse 
of a length). Combining Eqs. (|4T1) and (fTT]) yields 

uj ■ V0(r, v|ro, vo) + St(r, v)0(r, v|ro, vo) = 

= J dv'C(v|v', r)St(r, v')0(r, v'|ro, vq) + S (42) 

where ^'(r, v|ro, Vq) = i;t(r, v)0(r, v|ro, vq) [sS, [s^I- 
The quantity (/)(r, v|ro, vq) is called flux in neutronics. 

In the following, we introduce some simplifying hy- 
potheses. First, we consider a Id setup, where particles 
undergo exponential displacements along a straight line, 
only forward and backward directions being allowed: the 
so-called rod model [l], [3, El]. Further, we assume that 
the particle energy is preserved along each trajectory, 
which corresponds to the so-called one-speed approxima- 
tion. Finally, we take the total cross section and the scat- 
tering probability to be constant, i.e., 'Et{'r,v) = Ei and 
p(r,v) = p. Without loss of generality, we set St = 1. 
Despite being admittedly oversimplified, this model can 
nonetheless capture the essential features of the corre- 
sponding physical system. 

We define ujf and Wf, the forward and backward direc- 



tions, respectively. Similarly, we denote by 5/ and Sb the 
forward and backward intensities of the source, located 
at xq. Anisotropy is taken into account by introducing 
a persistence coefficient a such that after each scattering 
collision the particle preserves its direction uj with proba- 
bility a, and inverses its direction otherwise [H, |48[. This 
imposes 

Ciuj\uj') = p [a(5(w - uj') + (1 - a)6iuj + w')] • (43) 

The persistence coefficient is related to the average scat- 
tering angle by —1 -I- 2a = /z. Remark that a = 1 corre- 
sponds to a walker that systematically preserves its inci- 
dent direction (forward scattering), whereas a = to a 
walker that systematically reverses its incident direction. 
The model with a = has long been investigated under 
the name of telegrapher's equation 0, 111,143]. Isotropic 
scattering is recovered for a = 1/2 (/U = 0). For expo- 
nential flights 7 = 1/2, so that we have ij ~ 1 — fi. The 
volume V is assumed to be the interval V — [—R,R]. 
With this choice of parameters and notations, Eq. (j42|) 
yields the following set of stationary first-order differen- 
tial equations for the incident collision density 



J 



— + 1 j '^{x,ujf\xo,ujf) = p[a'^{x,ujf\xo,ujf) + (1 - a)'^{x,uJb\xo,ujf)] + Sf, 
+ 1^ ■i'{x,uJb\xo,ujf) = p [5'(x,Wb]xo,w/) + (1 - a)'if{x,ujf\xo,ujf)] , 



dx 



— + 1 j 'S{x,ujf\xo,ujb) = p [a'9{x,ujf\xo,ujb) + (1 - a)'<i'{x,uJb\xo,ujb)] , 
+ 1 I 'i'{x,ujb\xo,ujb) = p [a'i'{x,ujb\xo,ujb) + (1 - a)'i'{x,ujf\xo,uJb)] + Sb- 



(44) 



r 



These equations are linear and can be put in matricial 
form: together with appropriate boundary conditions, 
this leads to an explicit solution for 'i'{x,uj\xo, ujq). How- 
ever, the solution may not exist, depending on the choice 
of the equation parameters. Actually, Id exponential 
flights are recurrent wal ks ( i.e. . they almost surely re- 
visit their initial position [3i|,[3a]), so that "^{x,uj\xo,ujo) 
diverges when p = 1, unless leakage boundary conditions 
are imposed, i.e., the particles are lost upon crossing the 
boundary of the domain. 



B. Examples of calculations 



The effects of the scattering probabilit y p on the 



34] 



Here 



moments have been discussed elsewhere 
we will focus on the case of leakage boundary con- 
ditions without absorption, which allows emphasizing 



the effects of anisotropy. Leakages at a; = ±i? im- 
pose 'i'{-R,u!f\xo,ujf) = 0, *(i?, Wfc]a::o,a;/) ^ 0, 
"i'{R,uJb\xo,iUb) — 0, and '^{—R,ujf\xo,uJb) — 0, which 
corresponds to an homogeneous medium surrounded by 
vacuum. We choose to observe events when particles 
enter collisions, which amounts to counting scattering 
and absorptions, and disregarding the contribution of 
the source. Once '^(x,uj\xq,ujo) is known, the moments 
{ny)p{xo) can be obtained from Eq. (|^. For instance, 
the average ny reads 



(2:0) 



2R + T]R'^ 



rixl 



(45) 



when the source is isotropic. This expression explicitly 
depends on the initial position xq , on the persistence co- 
efficient a and on the size R of the domain, when = 1. 
Moments for generic St are simply obtained by rescaling 
the space variables R and xq by a factor Ef . In Fig. [T] 



7 



we display the behavior of (riy) p{xq) for R = 1. Remark 
that the average has always a maximum at = (whose 
height decreases with a) and a minimum at xq = R (in- 
dependent of a). When a = 1, the average collision num- 
ber becomes independent of the starting point xq (recall 
that we have chosen an isotropic source). The moment 
(ny)Q(xo) can be computed based on Eqs. and (ITOl) . 
and reads 



(ny)Q(xo) = {ny)p{xo) + 1. 



(46) 



When R is large (i.e., when i?Ei ^ 1), by comparing 
Eq. (j45|) with the results for isotropic transport in [3l| . 
we have {iiy) pixo)aniso — vinv) p{xo)iso, in agreement 
with Eq. (gni). 



The second moment can be also easily obtained by 
direct integration, and the formula reads 



(ny)p(xo) = {ny)p{xo) + 



24i?^ -I- 20r/i?^ -f bifR'^ - UrjRx^ 



(Srj^R^xl + ti'^Xq 



12 



+ {xl-R') 



(47) 



r 



In Fig. [2]we display the behavior of {ny)p{xo) ior R ~ I 
and an isotropic source. In this case, the curves are not 
monotonic with respect to xq: by varying a, the maxi- 
mum of (riy) p{xo) is either at the center of the domain, 
or at its boundary. The moment {ny)Q{xQ) can be com- 
puted based on Eqs. ([551) and (ITOl) . and reads 

{n^Qixa) = (n^)p(xo) + 2{nl,)p{xo) + 1. (48) 
Again, by comparison with the results for isotropic 



transport, when R is large we have 



{xo)a 



r]^{ny)p{xo)iso, in agreement with Eq. pO|) . 

Then, for the same geometrical configuration (i.e., 
leakage boundaries and R — 1), we analyze the behav- 
ior of the moments when observed up to entering the 
n-th collision, as a function of n for xq — and vary- 
ing a. These quantities, that we denote by {riy) p{n, xq), 
are easily obtained by Monte Carlo simulation, and are 
expected to converge to the respective (n™)p(xo) when 
n —>■ oo. Figures |3] and |4] show that the average and the 
second moment grow with n and saturate to their respec- 
tive asymptotic values, given by Eqs. ps)) and (|T7|) . The 
number of collisions needed to reach saturation, as well 
as the asymptotic value at saturation, decrease with in- 
creasing a. This can be understood by considering that 
a forward-peaked walker (a ~ 1) will undergo fewer col- 
lisions in V (before crossing the boundary) than a walker 
with small a, which on the contrary frequently reverses 
its direction and thus stays longer in V. 

Finally, we conclude with some considerations concern- 
ing the limit behavior of the collision number distribu- 
tion. To fix the ideas, let us assume that xq = 0. When 
R is large, by direct inspection we recognize the scaling 



'i?^™, where c„ 



(-l)"m!E2m/(2m)!, 



and 



E2r, 



2n+l k 
fc=l ]=Q 



-ly 



k\ [k - 2j) 



2ri+l 



k{2i) 



(49) 



function G{u) associated to P{nv) has the moment ex- 
pansion 



G{u) 



+ CXD 

E( 

m=0 



{-uf 



+ 00 

E (2m)\ 

m=0 ^ ' 



{R'ljur. (50) 



are the Euler's numbers [36|. The moment generating 



By carrying out the sum we get G{u) = sech{R^/rju), 
where we recognize the scaling variable yy = ny / iv^'^)- 
In Fig. [5] we display P{yv) as a function of yv, for var- 
ious values of rj and R: it is immediately apparent that 
all the curves collapse. The exponential tail of P{yv) 
for large W_is expected on the basis of the Tauberian 
theorems [5(l|, since G{u) ~ 1 — riR^u/2 for small u. 

All analytical results reported here have been validated 
by comparison with Monte Carlo simulations with 10^ 
particles. 



VI. CONCLUSIONS 

In this paper we have proposed a general approach to 
the counting statistics for the number of events falling 
within a given region V of the phase space, when the 
underlying process is a random flight. By resorting to 
survival probabilities, we have provided an explicit de- 
scription of the moments and factorial moments of the 
collision number ny . Only a minimal number of hypothe- 
ses on the underlying transport kernels K{z\z') or L(z|z') 
are required, and we were able to take into account the 
effects of jump length distribution, anisotropy, absorp- 
tion and boundary conditions. In this work we have in 
particular focused on anisotropy, and discussed some ex- 
amples of analytical calculations for the class of random 
flights where displacements are exponentially distributed, 
namely, the exponential flights. 

In view of the physical systems to which random flights 
most often apply, e.g., gas dynamics and neutronics, we 
have found natural here to resort to the language specific 
to stochastic transport phenomena. However, this same 
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formalism could be further generalized in terms of semi- 
Markov renewal processes for an arbitrary state variable 
q evolving in the phase space according to some transi- 
tion kernel. 

Finally, observe that we might take advantage of the 
knowledge on the number of collisions in a given domain, 
when available, as an estimator to infer the equilibrium 
distribution of the underlying stochastic path, which is 
seldom directly accessible. Indeed, while for instance in 
neutronics the underlying transport process is suppos- 
edly known, and one is typically interested in assessing 
the collision statistics (e.g., the deposited energy and/or 
the radiation damage) , in life sciences one could measure 
the number of hits of the displacing species in a region V 
so as to probe its possibly unknown dynamics. Nonethe- 
less, this inverse problem might be ill-posed, and deserves 
further investigation. 



Appendix A: Operator identities 

In this Section we recall some identities that are used 
in the paper. The formal Neumann series is defined by 



oo 

EA"[/](^) = f^[/](^): 



(Al) 



where I is the identity operator. 

For any sufficiently well-behaved operator A, we have 
the identity 



+00 



(l-A)^n™A"-M/](z) = 



n=0 



A 



A 



k-l 



im, (A2) 



where the coefficients 



i=0 



are the Stirling numbers of second kind 
we have 



(A3) 



Moreover, 



+ CX3 



I 



fc=0 n=0 



= m ' 



I- A 



im, (A4) 



where the coefficients Sm,fc are the Stirling numbers of 
the first kind [s^ , which are defined as the coefficients in 
the expansion 



(A5) 



fc=0 



{x)m being the falling factorial 

l){x ~ 2)...{x - m + 1). (A6) 
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